Violations of physical and psychological expectations in the human adult brain - Experiment 1 Univariate Analysis
Shari Liu
2023-08-01
1 Overview
This script:
- Part 1: visualizes and analyzes the mean amplitude of response
(
meanbeta) in 8 focal regions (left and right SMG, left and right STS, left and right V1, bilateral APC and RFC). - Part 2: does the same across other, non-focal ROIs.
- Part 3: visualizes and analyzes the reliability of domains with events and events within domains across all ROIs. This is the analog to the pre-registered MVPA analysis that is reported in the SI.
It also contains the code needed to reproduce Figure 3 in the main
text of the paper. All instances of write_rds,
write_csv, etc have been commented out to avoid writing
over existing files. Uncomment those statements in order to reproduce
these steps and save (new) outputs.
1.1 Read in data
univariate_data <- readRDS(here("outputs/univariate_data/study1_univariate_data_singlebetas.Rds"))
univariate_data$event <- relevel(as.factor(univariate_data$event), ref = "unexp")
univariate_summary_domain <- readRDS(here("outputs",
"univariate_descriptive_summaries",
"study1_univariate_summary_domain_allvoxelN.Rds"))
univariate_summary_task <- readRDS(here(
"outputs",
"univariate_descriptive_summaries",
"study1_univariate_summary_task_allvoxelN.Rds"))2 PART 1: Focal regions
focal_data <- univariate_data %>%
filter(focal_region == 1)
focal_summary_domain <- univariate_summary_domain %>%
filter(focal_region == 1)
focal_summary_task <- univariate_summary_task %>%
filter(focal_region == 1)
domain_specific_regions <- c("physics", "psychology")
domain_general_regions <- c("early_visual", "MD")focal_data_100 <- focal_data %>%
filter(top_n_voxels == "100")
focal_summary_domain_100 <- focal_summary_domain %>%
filter(top_n_voxels == "100")
focal_summary_task_100 <- focal_summary_task %>%
filter(top_n_voxels == "100")2.1 Does the VOE effect diminish in size over runs?
model_habituation <- lmer(data = focal_data_100 %>%
filter(ROI_category != "early_visual",
event != "fam",
ROI_category != "MD"), # cannot model MD regions because of non-independence of ROI selection (based on runs 2-4)
formula = meanbeta ~ extracted_run_number * event + (1|subjectID))
check_model(model_habituation)
plot(allEffects(model_habituation))
summary(model_habituation)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ extracted_run_number * event + (1 | subjectID)
## Data: focal_data_100 %>% filter(ROI_category != "early_visual", event !=
## "fam", ROI_category != "MD")
##
## REML criterion at convergence: 17809
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8142 -0.5931 -0.0313 0.5983 5.9668
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 0.2364 0.4862
## Residual 3.4448 1.8560
## Number of obs: 4352, groups: subjectID, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.181e-01 1.212e-01 1.600e+01 2.624 0.01843
## extracted_run_number1 5.679e-02 4.873e-02 4.328e+03 1.165 0.24396
## extracted_run_number2 -9.594e-02 4.873e-02 4.328e+03 -1.969 0.04904
## extracted_run_number3 7.969e-03 4.873e-02 4.328e+03 0.164 0.87010
## event1 7.583e-02 2.814e-02 4.328e+03 2.695 0.00706
## extracted_run_number1:event1 1.519e-01 4.873e-02 4.328e+03 3.117 0.00184
## extracted_run_number2:event1 -1.063e-02 4.873e-02 4.328e+03 -0.218 0.82738
## extracted_run_number3:event1 -6.623e-02 4.873e-02 4.328e+03 -1.359 0.17417
##
## (Intercept) *
## extracted_run_number1
## extracted_run_number2 *
## extracted_run_number3
## event1 **
## extracted_run_number1:event1 **
## extracted_run_number2:event1
## extracted_run_number3:event1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ext__1 ext__2 ext__3 event1 e__1:1 e__2:1
## extrctd_r_1 0.000
## extrctd_r_2 0.000 -0.333
## extrctd_r_3 0.000 -0.333 -0.333
## event1 0.000 0.000 0.000 0.000
## extrct__1:1 0.000 0.000 0.000 0.000 0.000
## extrct__2:1 0.000 0.000 0.000 0.000 0.000 -0.333
## extrct__3:1 0.000 0.000 0.000 0.000 0.000 -0.333 -0.333model_habituation_results <- cbind(gen.m(model_habituation), gen.ci(model_habituation)[3:10,])## Computing profile confidence intervals ...model_habituation_results_byrun <- lsmeans(model_habituation, pairwise ~ event | extracted_run_number, pbkrtest.limit = 4352)$contrasts %>% as.data.frame()DT::datatable(model_habituation_results, options = list(scrollX = TRUE))DT::datatable(model_habituation_results_byrun, options = list(scrollX = TRUE))# write_rds(model_habituation_results, path = here("outputs/univariate_results/habituation_over_runs_modelsummary.Rds"))
# write_rds(model_habituation_results_byrun, path = here("outputs/univariate_results/habituation_over_runs_perrun.Rds"))2.2 Visualize by domain
regions <- levels(as.factor(focal_data$ROI_name_final))
plot_by_domain_run1 <- function(region, ymin, ymax) {
plotobject <- ggplot(data = focal_summary_domain_100 %>% filter(str_detect(ROI_name_final, region), extracted_run_number == "run1", event != "fam"),
aes(x = event, y = meanbeta, fill = domain)) +
geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
geom_errorbar(
aes(ymin = meanbeta - se, ymax = meanbeta + se),
position = position_dodge(width = .9),
width = .2,
colour = "black"
) +
theme_cowplot(10) +
facet_wrap( ~ ROI_name_final + domain, nrow = 1) +
scale_fill_manual(values = c("#00AFBB", "#FC4E07")) +
ylab("Amplitude") +
xlab("Event") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
),
legend.position = "none") +
coord_cartesian(ylim = c(ymin, ymax)) +
ggtitle(paste0("ROI:", region))
plotobject
}This plot is a subset of Figure 3 in the paper.
SMG <- plot_by_domain_run1("supramarginal_", -1, 4)
STS <- plot_by_domain_run1("superiortemporal_", 0, 4)
early_vis <- plot_by_domain_run1("V1_", 0, 11)
MD_APC <- plot_by_domain_run1("antParietal_bilateral", -1, 1)
MD_RFC <- plot_by_domain_run1("precentral_inferiorfrontal_R", -1, 1)
(SMG + STS) / (early_vis + MD_APC + MD_RFC) + plot_layout(widths = c(1, 1, 2, 1))## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
2.2.1 Plot run by run (except for MD regions)
regions <- levels(as.factor(focal_data$ROI_name_final))
plot_by_domain_all_runs <- function(region, ymin, ymax) {
plotobject <- ggplot(data = focal_summary_domain_100 %>% filter(str_detect(ROI_name_final, region)),
aes(x = event, y = meanbeta, fill = domain)) +
geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
geom_errorbar(
aes(ymin = meanbeta - se, ymax = meanbeta + se),
position = position_dodge(width = .9),
width = .2,
colour = "black"
) +
theme_cowplot(10) +
facet_wrap( ~ ROI_name_final + top_n_voxels + domain + extracted_run_number, nrow = 1) +
scale_fill_manual(values = c("#00AFBB", "#FC4E07")) +
ylab("Average beta (amplitude)") +
xlab("Event") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
),
legend.position = "none") +
coord_cartesian(ylim = c(ymin, ymax)) +
ggtitle(paste0("ROI:", region))
plotobject
}plot_by_domain_all_runs("superiortemporal_L", -1, 3) | plot_by_domain_all_runs("superiortemporal_R", -1, 3)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_by_domain_all_runs("supramarginal_L", -1, 1) | plot_by_domain_all_runs("supramarginal_R", -1, 1)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_by_domain_all_runs("V1_L", -1, 7) | plot_by_domain_all_runs("V1_R", -1, 7)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_by_domain_run1("antParietal_bilateral", -1, 1) | plot_by_domain_run1("precentral_inferiorfrontal_R", -1, 1)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
2.3 Visualize by task
regions <- levels(as.factor(focal_summary_task_100$ROI_name_final))
plot_univar_bytask <- function(region, ymin, ymax) {
plotobject <- ggplot(data = focal_summary_task_100 %>%
filter(str_detect(ROI_name_final, region),
extracted_run_number == "run1"),
aes(x = event, y = meanbeta, fill = task)) +
geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
geom_errorbar(
aes(ymin = meanbeta - se, ymax = meanbeta + se),
position = position_dodge(width = .9),
width = .2,
colour = "black"
) +
# geom_point(data = focal_data_100 %>% filter(str_detect(ROI_name_final, region),
# extracted_run_number == "run1"),
# alpha = .1) +
# geom_line(data = focal_data_100 %>% filter(str_detect(ROI_name_final, region)),
# aes(group =
# subjectID),
# alpha = .1) +
theme_cowplot(10) +
facet_wrap( ~ ROI_name_final + factor(task, c("solidity", "support", "goal", "efficiency")), nrow = 1) +
scale_fill_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
ylab("Average beta (amplitude)") +
xlab("Event") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
),
legend.position = "none") +
coord_cartesian(ylim = c(ymin, ymax)) +
ggtitle(paste0("ROI:", region))
plotobject
}plot_univar_bytask("superiortemporal_L", -1, 3) | plot_univar_bytask("superiortemporal_R", -1, 3)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_bytask("supramarginal_L", -1, 3) | plot_univar_bytask("supramarginal_R", -1, 3)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_bytask("V1_L", -1, 7) | plot_univar_bytask("V1_R", -1, 7)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
plot_univar_bytask("antParietal_bilateral", -1, 3) | plot_univar_bytask("precentral_inferiorfrontal_R", -1, 3)## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
2.4 Results
2.4.1 Model fitting per domain
effect_sizes_intercept_row <- data.frame(d = NA)
BF_intercept_row <- data.frame(BF = NA, BF_interpretation = NA)
rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames(BF_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1", "domain1", "rep.L", "event:domain")
data_formodeling <- univariate_data %>%
filter(event != "fam")
data_formodeling %>%
group_by(top_n_voxels) %>%
summarise(n = n())
ROIs <- levels(as.factor(data_formodeling$ROI_name_final))
top_n_voxels <- unique(data_formodeling$top_n_voxels)
runs <- c("run1", "run2", "run3", "run4", "all_runs")
modelsummaries_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
for (topN in top_n_voxels) {
cur_top_voxels_N <- topN
for (run in runs) {
curRun <- run
if (curRun != "all_runs") {
ROIdata <- data_formodeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN,
extracted_run_number == curRun)
model <- lmer(data = ROIdata,
formula = meanbeta ~ event * domain + rep + (1 |
subjectID))
} else if (curRun == "all_runs") {
ROIdata <- data_formodeling %>%
filter(ROI_name_final == curROI,
top_n_voxels == topN)
model <- lmer(data = ROIdata,
formula = meanbeta ~ event * domain + rep + (1 |
subjectID))
}
model_just_maineffects <-
update(model, formula = ~ . - event:domain) # Without interaction term
model_just_domain <-
update(model_just_maineffects, formula = ~ . - event)
model_just_event <-
update(model_just_maineffects, formula = ~ . - domain)
model_no_rep <- update(model, formula = ~ . - rep)
BF_BIC_interaction <-
data.frame(BF = exp((
BIC(model_just_maineffects) - BIC(model)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_maineffects) - BIC(model)
) / 2))) # BICs to Bayes factor
BF_BIC_event <-
data.frame(BF = exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_domain) - BIC(model_just_maineffects)
) / 2)))
BF_BIC_domain <-
data.frame(BF = exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_just_event) - BIC(model_just_maineffects)
) / 2)))
BF_BIC_rep <-
data.frame(BF = exp((BIC(model_no_rep) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_no_rep) - BIC(model)
) / 2)))
BF_initial <- rbind(BF_BIC_event,
BF_BIC_domain,
BF_BIC_rep,
BF_BIC_interaction)
rownames(BF_initial) <- rownames_no_intercept
BFs <- rbind(BF_intercept_row,
BF_initial)
effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
select(d)
effect_sizes <-
rbind(effect_sizes_intercept_row, effect_sizes_initial)
if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
ROI_category <- "domain_specific"
} else {
ROI_category <- "domain_general"
}
ROI_focal <- ROIdata$focal_region[1]
modelsummary <-
cbind(
cur_top_voxels_N,
curRun,
ROI_category,
ROI_focal,
curROI,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:7, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"d",
"BF",
"BF_interpretation"
)
modelsummaries_init <-
rbind(modelsummaries_init, modelsummary)
cat(
tab_model(
model,
show.std = TRUE,
show.se = TRUE,
show.stat = TRUE,
auto.label = TRUE,
title = paste0("Univariate ROI results in: ", curROI)
)$knitr,
" \n \n"
)
}
}
}modelsummaries <- modelsummaries_init %>%
mutate_at(c(7:15), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
modelsummaries$effect <-
as.factor(modelsummaries$effect)
levels(modelsummaries$effect) <-
c("intercept", "domain", "event", "event:domain", "rep")
rownames(modelsummaries) <- NULL# write_rds(modelsummaries, here("outputs", "univariate_results", "modelsummaries_allregions_alltopN.Rds"))focal_modelsummaries <- read_rds(here("outputs/univariate_results/modelsummaries_allregions_alltopN.Rds")) %>%
filter(focal_region == 1)
DT::datatable(focal_modelsummaries %>% filter(top_n_voxels == "100", focal_region == 1) %>% arrange(effect), options = list(scrollX = TRUE,
pageLength = 100)) # filter(extracted_run_number == "all_runs" | extracted_run_number == "run1") %>%2.4.1.1 LSMG
LSMG_model <- lmer(data = univariate_data %>%
filter(top_n_voxels == "100",
event != "fam",
extracted_run_number == "run1",
ROI_name_final == "supramarginal_L"),
formula = meanbeta ~ event * domain + rep + (1|subjectID))
LSMG_results_bydomain_top100 <- lsmeans(LSMG_model, pairwise ~ event | domain)$contrasts %>% as.data.frame()
DT::datatable(LSMG_results_bydomain_top100, options = list(scrollX = TRUE))# write_rds(LSMG_results_bydomain_top100, path = here("outputs/univariate_results/LSMG_results_bydomain_top100.Rds"))2.4.1.2 RSMG
RSMG_model <- lmer(data = univariate_data %>%
filter(top_n_voxels == "100",
event != "fam",
extracted_run_number == "run1",
ROI_name_final == "supramarginal_R"),
formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(RSMG_model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
## Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",
## extracted_run_number == "run1", ROI_name_final == "supramarginal_R")
##
## REML criterion at convergence: 1047.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2569 -0.5727 0.0203 0.6278 2.5903
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.517 1.232
## Residual 2.312 1.520
## Number of obs: 272, groups: subjectID, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.01345 0.31266 16.00000 -0.043 0.96621
## event1 0.21889 0.09219 251.00000 2.374 0.01834 *
## domain1 0.28691 0.09219 251.00000 3.112 0.00207 **
## rep.L -0.36357 0.13038 251.00000 -2.789 0.00570 **
## event1:domain1 0.05727 0.09219 251.00000 0.621 0.53505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 rep.L
## event1 0.000
## domain1 0.000 0.000
## rep.L 0.000 0.000 0.000
## event1:dmn1 0.000 0.000 0.000 0.000lsmeans(RSMG_model, pairwise ~ event | domain)$contrasts## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.552 0.261 251 2.118 0.0352
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.323 0.261 251 1.240 0.2163
##
## Results are averaged over the levels of: rep
## Degrees-of-freedom method: kenward-roger2.4.1.3 LSTS
LSTS_model <- lmer(data = univariate_data %>%
filter(top_n_voxels == "100",
event != "fam",
extracted_run_number == "run1",
ROI_name_final == "superiortemporal_L"),
formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(LSTS_model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
## Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",
## extracted_run_number == "run1", ROI_name_final == "superiortemporal_L")
##
## REML criterion at convergence: 1051.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4985 -0.5627 -0.0406 0.6461 2.9036
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 1.773 1.331
## Residual 2.327 1.525
## Number of obs: 272, groups: subjectID, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.54190 0.33591 16.00000 1.613 0.126244
## event1 0.16841 0.09249 251.00000 1.821 0.069839 .
## domain1 -0.32878 0.09249 251.00000 -3.555 0.000452 ***
## rep.L 0.13244 0.13081 251.00000 1.012 0.312286
## event1:domain1 0.04951 0.09249 251.00000 0.535 0.592897
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 rep.L
## event1 0.000
## domain1 0.000 0.000
## rep.L 0.000 0.000 0.000
## event1:dmn1 0.000 0.000 0.000 0.000lsmeans(LSTS_model, pairwise ~ event | domain)$contrasts## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.436 0.262 251 1.666 0.0970
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.238 0.262 251 0.909 0.3643
##
## Results are averaged over the levels of: rep
## Degrees-of-freedom method: kenward-roger2.4.1.4 RSTS
RSTS_model <- lmer(data = univariate_data %>%
filter(top_n_voxels == "100",
event != "fam",
extracted_run_number == "run1",
ROI_name_final == "superiortemporal_R"),
formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(RSTS_model)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
## Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",
## extracted_run_number == "run1", ROI_name_final == "superiortemporal_R")
##
## REML criterion at convergence: 1015.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.16150 -0.59428 -0.03496 0.59256 2.87577
##
## Random effects:
## Groups Name Variance Std.Dev.
## subjectID (Intercept) 2.143 1.464
## Residual 1.992 1.411
## Number of obs: 272, groups: subjectID, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.95600 0.36522 16.00000 2.618 0.018662 *
## event1 0.25840 0.08557 251.00000 3.020 0.002792 **
## domain1 -0.30885 0.08557 251.00000 -3.609 0.000371 ***
## rep.L -0.22514 0.12102 251.00000 -1.860 0.064000 .
## event1:domain1 -0.05083 0.08557 251.00000 -0.594 0.553030
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) event1 doman1 rep.L
## event1 0.000
## domain1 0.000 0.000
## rep.L 0.000 0.000 0.000
## event1:dmn1 0.000 0.000 0.000 0.000lsmeans(RSTS_model, pairwise ~ event | domain)$contrasts## domain = physics:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.415 0.242 251 1.715 0.0876
##
## domain = psychology:
## contrast estimate SE df t.ratio p.value
## unexp - exp 0.618 0.242 251 2.555 0.0112
##
## Results are averaged over the levels of: rep
## Degrees-of-freedom method: kenward-roger2.4.2 Model fitting per task
effect_sizes_intercept_row <- data.frame(d = NA)
BF_intercept_row <- data.frame(BF = NA, BF_interpretation = NA)
rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames(BF_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1", "rep.L")
data_formodeling <- univariate_data %>%
filter(event != "fam",
top_n_voxels == 100)
data_formodeling %>%
group_by(top_n_voxels) %>%
summarise(n = n())
ROIs <- levels(as.factor(data_formodeling$ROI_name_final))
top_n_voxels <- unique(data_formodeling$top_n_voxels)
runs <- c("run1", "all_runs")
tasks <- unique(data_formodeling$task)
ROIs <- unique(data_formodeling$ROI_name_final)
data_formodeling$event <- relevel(data_formodeling$event, ref = "unexp")
modelsummaries_pertask_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
for (task in tasks) {
cur_task <- task
for (run in runs) {
if (run == "run1") {
ROIdata <-
data_formodeling %>% filter(ROI_name_final == curROI,
task == cur_task,
extracted_run_number == "run1")
} else if (run == "all_runs")
ROIdata <-
data_formodeling %>% filter(ROI_name_final == curROI,
task == cur_task)
model <- lmer(data = ROIdata,
formula = meanbeta ~ event + rep + (1 | subjectID)) # had to remove run random int (convergence))
focal_region <- ROIdata$focal_region[1]
task_domain <- as.character(ROIdata$domain[1])
ROI_category <- ROIdata$ROI_category[1]
model_null <-
update(model, formula = ~ . - event) # Without interaction term
model_norep <-
update(model, formula = ~ . - rep) # Without interaction term
BF_BIC_event <-
data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_null) - BIC(model)
) / 2))) # BICs to Bayes factor
BF_BIC_rep <-
data.frame(BF = exp((BIC(model_norep) - BIC(model)) / 2),
BF_interpretation = interpret_bf(exp((
BIC(model_norep) - BIC(model)
) / 2))) # BICs to Bayes factor
# BF_initial <- rbind(BF_BIC_event)
BFs <- rbind(BF_intercept_row,
BF_BIC_event,
BF_BIC_rep)
rownames(BFs)[2:3] <- rownames_no_intercept
effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
select(d)
effect_sizes <-
rbind(effect_sizes_intercept_row, effect_sizes_initial)
cur_top_voxels_N <- 100
modelsummary <-
cbind(
cur_top_voxels_N,
run,
ROI_category,
focal_region,
curROI,
cur_task,
rownames(summary(model)$coefficients),
summary(model)$coefficients,
confint.merMod(model)[3:5, ],
effect_sizes,
BFs
) %>%
as.data.frame()
colnames(modelsummary) <-
c(
"top_n_voxels",
"extracted_run_number",
"domain",
"focal_region",
"ROI",
"task",
"effect",
"B",
"SE",
"df",
"t",
"p",
"CI_95_lower",
"CI_95_upper",
"d",
"BF",
"BF_interpretation"
)
modelsummaries_pertask_init <-
rbind(modelsummaries_pertask_init, modelsummary)
# model_parameters(model) %>% print_md()
cat(
tab_model(
model,
show.std = TRUE,
show.se = TRUE,
show.stat = TRUE,
auto.label = TRUE,
title = paste0(
"Exploratory Univariate ROI results for the task ",
cur_task,
"in the region: ",
curROI,
"in ",
curRun
)
)$knitr,
" \n \n"
)
}
}
}modelsummaries_pertask <- modelsummaries_pertask_init %>%
mutate_at(c(8:16), as.numeric) %>%
mutate(star = as.factor(
case_when(p < .001 ~ "***",
p < .01 ~ "**",
p < .05 ~ "*",
p < .1 ~ "~",
TRUE ~ " ")
)) %>%
mutate(p = round(p, digits = 3))
modelsummaries_pertask$effect <-
as.factor(modelsummaries_pertask$effect)
levels(modelsummaries_pertask$effect) <-
c("intercept", "event", "rep")
rownames(modelsummaries_pertask) <- NULL# write_rds(modelsummaries_pertask, here("outputs", "univariate_results", "modelsummaries_allregions_pertask.Rds"))modelsummaries_pertask <- read_rds(here("outputs", "univariate_results", "modelsummaries_allregions_pertask.Rds"))2.4.2.1 Plot model outputs
focal_modelsummaries_pertask_forplotting <- modelsummaries_pertask %>%
filter(top_n_voxels == "100",
effect == "event",
extracted_run_number == "run1",
focal_region == 1)
focal_modelsummaries_pertask_forplotting$ROI <- factor(focal_modelsummaries_pertask_forplotting$ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "antParietal_bilateral", "precentral_inferiorfrontal_R", "V1_L", "V1_R"))
pertask_b <- ggplot(focal_modelsummaries_pertask_forplotting,
aes(ROI, B, colour = task, group = task)) +
geom_point() +
geom_errorbar(
aes(ymin = B - SE, ymax = B + SE,
colour = task),
position = position_dodge(width = .9),
width = .2
) +
# geom_line() +
geom_hline(yintercept = 0) +
facet_wrap( ~ task, nrow = 1) +
scale_colour_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
ylab("B and SE") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) +
coord_cartesian(ylim = c(-.5, .7))
perregion_b <- ggplot(focal_modelsummaries_pertask_forplotting,
aes(task, B, colour = task, group = task)) +
geom_point() +
geom_errorbar(
aes(ymin = B - SE, ymax = B + SE,
colour = task),
position = position_dodge(width = .9),
width = .2
) +
# geom_line() +
geom_hline(yintercept = 0) +
facet_wrap( ~ ROI, nrow = 1) +
scale_colour_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
ylab("B and SE") +
theme(axis.text.x = element_text(
angle = 90,
vjust = 0.5,
hjust = 1
)) +
coord_cartesian(ylim = c(-.5, .7))pertask_b
perregion_b
2.5 Sensitivity of focal region results to top voxel criteria
focal_modelsummaries$top_n_voxels <- factor(focal_modelsummaries$top_n_voxels, levels = c("10", "50", "100", "150", "200", "250", "300"))
focal_modelsummaries$extracted_run_number <- factor(focal_modelsummaries$extracted_run_number, levels = c("run1", "run2", "run3", "run4", "all_runs"))
ggplot(focal_modelsummaries %>%
filter(extracted_run_number != "all_runs",
ROI != "antParietal_bilateral",
ROI != "precentral_inferiorfrontal_R",
effect != "intercept",
effect != "rep"),
aes(extracted_run_number, d, colour = top_n_voxels)) +
geom_point() +
geom_line(aes(group = top_n_voxels)) +
ylab("Cohen's D") +
facet_wrap(~ effect + factor(ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "V1_L", "V1_R")), nrow = 3) +
# coord_cartesian(ylim = c(-.75, 0.75)) +
geom_hline(yintercept = 0) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
3 PART 2: Other regions
nonfocal_modelsummaries <- read_rds(here("outputs/univariate_results/modelsummaries_allregions_alltopN.Rds")) %>%
filter(focal_region == 0)
DT::datatable(nonfocal_modelsummaries %>%
filter(top_n_voxels == "100") %>%
arrange(effect),
options = list(scrollX = TRUE,
pageLength = 100))4 PART 3: Region by region univar effect size
Here we repeat the MVPA manyregions analysis, whose results are difficult to interpret given the lack of reliable multivariate signal towards unexpected events, despite clear univariate effects.
data_for_modeling <- univariate_data %>%
filter(event != "fam",
top_n_voxels == "100") %>%
filter(extracted_run_number == "run1") %>%
filter(!str_detect(ROI_name_final, "bilateral"))
ROIs <- levels(as.factor(data_for_modeling$ROI_name_final))
effect_sizes_intercept_row <- data.frame(d = NA)
rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1")modelsummaries_init <- data.frame()
for (ROI in ROIs) {
curROI <- ROI
# FIRST COMPUTE EVENT EFFECT SIZES, FOR EACH DOMAIN
physics_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
domain == "physics")
psychology_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
domain == "psychology")
# had to remove random intercept for runs (convergence))
physics_model <- lmer(data = physics_data,
formula = meanbeta ~ event + iqr_contrast + mean_luminance + mean_motion + mean_hsf + mean_lsf + mean_rect + mean_curv + (1|subjectID))
psychology_model <- lmer(data = psychology_data,
formula = meanbeta ~ event + (1|subjectID))
physics_event_effect_size <- lme.dscore(physics_model, physics_data, "lme4") %>%
select(d)
psychology_event_effect_size <- lme.dscore(psychology_model, psychology_data, "lme4") %>%
select(d)
# SECOND COMPUTE DOMAIN EFFECT SIZES, FOR EACH EVENT
expected_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
event == "exp")
unexpected_data <- data_for_modeling %>%
filter(ROI_name_final == curROI,
event == "unexp")
# had to remove random intercept for runs (convergence))
exp_model <- lmer(data = expected_data,
formula = meanbeta ~ domain + (1 | subjectID))
unexp_model <- lmer(data = unexpected_data,
formula = meanbeta ~ domain + (1 | subjectID))
exp_domain_effect_size <-
lme.dscore(exp_model, expected_data, "lme4") %>%
select(d)
unexp_domain_effect_size <-
lme.dscore(unexp_model, unexpected_data, "lme4") %>%
select(d)
if (physics_data$ROI_category[1] == "MD" | physics_data$ROI_category[1] == "early_visual") {
ROI_domain <- "general"
} else {
ROI_domain <- "specific"
}
modelsummary <-
cbind(
ROI_domain,
physics_data$ROI_category[1],
curROI,
physics_event_effect_size,
psychology_event_effect_size,
exp_domain_effect_size,
unexp_domain_effect_size
) %>%
as.data.frame()
modelsummaries_init <-
rbind(modelsummaries_init, modelsummary)
}colnames(modelsummaries_init) <-
c("ROI_domain",
"ROI_category",
"ROI",
"physics_event_effect_size",
"psychology_event_effect_size",
"exp_domain_effect_size",
"unexp_domain_effect_size"
)
rownames(modelsummaries_init) <- NULL
# write_rds(modelsummaries_init, here("outputs/univariate_results/univar_manyregions_results.Rds"))univar_manyregions_results <- read_rds( here("outputs/univariate_results/univar_manyregions_results.Rds"))plot1 <- ggplot(univar_manyregions_results, aes(psychology_event_effect_size, physics_event_effect_size)) +
geom_hline(yintercept = 0, linetype = 'dotted') +
geom_vline(xintercept = 0, linetype = 'dotted') +
coord_cartesian(xlim = c(-0.5, 1), ylim = c(-0.5, 1)) +
geom_point(aes(colour = ROI_category), size = 3) +
geom_smooth(method = "lm") +
facet_wrap(~ROI_domain, scales = "free") +
# coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
xlab("VOE Effect Size (d), Physics Events") +
ylab("VOE Effect Size (d), Psychology Events") +
ggtitle(label = "Exp 1 VOE effects for each domain across regions") +
geom_text_repel(aes(label = ROI), size = 3) +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.position = "none")
plot2 <- ggplot(univar_manyregions_results, aes(exp_domain_effect_size, unexp_domain_effect_size)) +
geom_hline(yintercept = 0, linetype = 'dotted') +
geom_vline(xintercept = 0, linetype = 'dotted') +
coord_cartesian(xlim = c(-1.2, 1.2), ylim = c(-1.2, 1.2)) +
geom_point(aes(colour = ROI_category), size = 3) +
geom_smooth(method = "lm") +
facet_wrap(~ROI_domain, scales = "free") +
# coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
xlab("Domain Effect Size (d), Expected Events") +
ylab("Domain Effect Size (d), Unexpected Events") +
ggtitle(label = "Exp 1 Domain effects for each event across regions") +
geom_text_repel(aes(label = ROI), size = 3) +
theme(axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
legend.position = "none")plot1## `geom_smooth()` using formula = 'y ~ x'
plot2## `geom_smooth()` using formula = 'y ~ x'
# instead of testing whether the linear relationship between x and y is 0,
# we test for independence instead. H0 is that x and y are independent; F_{XY}(x,y) = F_X(x) F_Y(y).
observed_cors_ind <- univar_manyregions_results %>%
group_by(ROI_domain) %>%
summarise(
cor_domain_within_event =
np.cor.test(
exp_domain_effect_size,
unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE
)$estimate,
p_domain_within_event =
np.cor.test(
exp_domain_effect_size,
unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE
)$p.value,
cor_event_within_domain =
np.cor.test(
physics_event_effect_size,
psychology_event_effect_size,
alternative = "two.sided",
independent = TRUE
)$estimate,
p_event_within_domain =
np.cor.test(
physics_event_effect_size,
psychology_event_effect_size,
alternative = "two.sided",
independent = TRUE
)$p.value
)observed_cors_ind# write_rds(observed_cors_ind, path = here("outputs/univariate_results/study1_univar_manyregions_observed_cor.Rds"))4.0.1 Bootstrap stats
set.seed(2020)
## First define a function to work out the difference of corrs:
diff_corr <- function(data, indices) {
data <- data[indices,]
cor1 <-
np.cor.test(DS_results$exp_domain_effect_size,
DS_results$unexp_domain_effect_size,
alternative = "two.sided",
independent = TRUE,
parallel = TRUE)$estimate
cor2 <-
np.cor.test(
data$psychology_event_effect_size,
data$physics_event_effect_size,
alternative = "two.sided",
independent = TRUE,
parallel = TRUE)$estimate
return(cor1 - cor2)
}DS_results <- univar_manyregions_results %>%
filter(ROI_domain == "specific")
DG_results <- univar_manyregions_results %>%
filter(ROI_domain == "general")4.0.2 Domain specific
doMC::registerDoMC(cores = 2) # for running in parallel
# Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DS <- boot(
# data = DS_results,
# R = 4000,
# statistic = diff_corr,
# stype = "i"
# )
# saveRDS(res_boot_DS, here("outputs/univariate_results/study1_univariate_dsregions_4000perms_confirmatory.Rds"))
res_boot_DS <- read_rds(here("outputs/univariate_results/study1_univariate_dsregions_4000perms_confirmatory.Rds"))
## Retrieve the empirical 95% confidence interval:
ds_results <- append(boot.ci(res_boot_DS, type = "perc", conf = 0.95),boot.pval(res_boot_DS))
# saveRDS(ds_results, here("outputs/univariate_results/study1_univariate_dsregions_summary.Rds"))4.0.3 Domain general
## Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DG <- boot(data = DG_results,
# R = 4000,
# statistic = diff_corr,
# stype = "i")
# saveRDS(res_boot_DG, here("outputs/univariate_results/study1_univariate_dgregions_4000perms_confirmatory.Rds"))
res_boot_DG <- readRDS(here("outputs/univariate_results/study1_univariate_dgregions_4000perms_confirmatory.Rds"))
plot(res_boot_DG)
## Retrieve the empirical 95% confidence interval:
dg_results <- append(boot.ci(res_boot_DG, type = "perc", conf = 0.95),boot.pval(res_boot_DG))
# saveRDS(dg_results, here("outputs/univariate_results/study1_univariate_dgregions_summary.Rds"))